Replica exchange molecular dynamics simulations of amyloid peptide aggregation 
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The replica exchange molecular dynamics (REMD) approach is applied to four oligomeric peptide 
systems. At physiologically relevant temperature values REMD samples conformation space and 
aggregation transitions more efficiently than constant temperature molecular dynamics (CTMD). 
During the aggregation process the energetic and structural properties are essentially the same 
in REMD and CTMD. A condensation stage toward disordered aggregates precedes the /3-sheet 
formation. Two order parameters, borrowed from anisotropic fluid analysis, are used to monitor the 
aggregation process. The order parameters do not depend on the peptide sequence and length and 
therefore allow to compare the amyloidogenic propensity of different peptides. 

Keywords: molecular dynamics; replica exchange; peptide aggregation; amyloid; order transition; nematic 
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INTRODUCTION 

A thorough samphng of conformational space is re- 
quired to describe the thermodynamics of complex sys- 
tems such as multiple peptide chains at finite con- 
centrations. Constant temperature molecular dynam- 
ics (CTMD) techniques often fail to adequately sample 
conformational space of frustrated and minimally frus- 
trated systems which are characterized by a rugged free- 
energy landscape where energy barriers between min- 
ima are higher than the thermal energy at physiolog- 
ical temperature. For this reason, a number of ap- 
proaches to enhance sampling of phase space have been 
introduced 0, IS B- The parallel tempering tech- 
nique (also known as replica exchange) was developed 
for dealing with the slow dynamics of disordered spin sys- 
tems 0. Sugita and Okamoto have extended the orig- 
inal formulation of replica exchange into an MD based 
version (REMD) and tested it on the pentapeptide Met- 
enkephalin in vacuo |^ . Although in the context of frag- 
ile liquids De Michele and Sciortino found that parallel 
tempering does not increase the speed of equilibration 
of the (slow) configurational degrees of freedom 0, in 
the case of atomistic simulations of proteins many differ- 
ent applications have shown the efficiency of the method. 
Sanbonmatsu and Garcia have used REMD to investi- 
gate the structure of Met-enkephalin in explicit water Q , 
and the a-helical stabilization by the arginine side-chain 
which was found to originate from the shielding of main 
chain hydrogen bonds [9| . REMD has also been applied 
to investigate the energy landscape of the C-terminal (3- 
hairpin of protein G [3 ^| and a three-helix bundle 
protein REMD in implicit solvent has been used to 
investigate the thermodynamics of designed 20-residue 
structured peptides (IJ, ll^ > and recently to study fold- 
ing of a helical transmembrane protein [l5| . 

Highly ordered protein aggregates are associated with 
severe human disorders including Alzheimer's disease, 



type II diabetes, systemic amyloidosis, and transmissi- 
ble spongiform encephalopathies Q. The soluble 
precursors of the ordered protein deposits do not share 
any sequence homology or common fold. However, X-ray 
diffraction data indicate a cross-/3-structure for most fib- 
rillar aggregates These findings suggest that key 
steps in the aggregation process may be common to all 
amyloidogenic proteins. Despite the medical relevance 
of amyloidoses, many important questions about the for- 
mation of ordered aggregates remain unanswered. There 
is experimental evidence that cytotoxicity is more pro- 
nounced for the early ag greg ates than for highly orga- 
nized fibrillar structures [2(]j. Moreover, some peptide 
fragments of amyloidogenic proteins display the same 
properties as the full-length protein, including cooper- 
ative kinetics of aggregation, fibril formation, binding 
of the dye Congo red, and the cross-/3 X-ray diffraction 
pattern |2lj |. Both findings are particularly interesting 
because current simulation approaches allow significant 
sampling only for oligomeric peptide systems. 

There have been several lattice studies on aggrega- 
tion in proteins. These simplified models have allowed 
to investigate the foldability and aggregation propen- 
sity 123, ii23j and how interaction potentials affect the 
properties of aggregation-prone proteins 24]. Harrison 
et al. have shown that less stable proteins have a greater 
chance of assuming alternative native states as multi- 
mers . MD simulations of aggregation have been per- 
formed by using a three-bead backbone and single-bead 
side chain model [2^ . While this simplified model has al- 
lowed the simulation of the competition between folding 
and aggregation for two four-helix bundles it is probably 
not possible to extract detailed information on energet- 
ics and sequence dependence. Recently, a minimalist Go 
model of four peptide strands has been investigated 
by MD simulations in a confining sphere and the aggre- 
gation process was shown to depend on both sequence 
and environment |28l |. Atomic models of amyloidogenic 
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peptides have been simulated 
treatment of the solvent [2 
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Recently, a replica exchange Monte Carlo technique 
has been applied to a lattice Go model of a minimalist 
multichain system to study the interplay between fold- 
ing and disordered aggregation but atomic model 
REMD applications to ordered aggregation have not been 
reported yet. 

In the present paper, REMD with implicit solvent 
is used to investigate the thermodynamics of the early 
steps of peptide aggregation and comparison is made with 
CTMD. The present work was motivated by three ques- 
tions: Is it possible to sample the early events of ordered 
peptide aggregation at physiologically relevant tempera- 
tures? Do the aggregation energetics sampled by REMD 
correspond to those observed in CTMD simulations? Are 
the nematic and polar order parameters, borrowed from 
liquid crystal theory, useful to describe aggregation? The 
simulation results indicate that all questions can be an- 
swered affirmatively. Moreover, the "liquid crystals" or- 
der parameters allow to discriminate amyloidogenic pep- 
tide sequences from those that form only disordered ag- 
gregates. 



for a molecular system with atomic nuclei located at r : 
(ri, ...jTat). The PARAM19 vacuo energy function is 
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where 6 is a bond length, 9 a bond angle, (p a dihedral 
angle, w an improper dihedral, r.^ is the distance between 
atoms i and j, qi and qj are partial charges, and d™'" and 
e™™ are the optimal van der Waals distance and energy, 
respectively. Parameters are given in Ref. [s^ . 

An implicit model based on the solvent accessible sur- 
face was used to describe the main effects of the aqueous 
solvent on the solute [s^l- In this approximation, the 
solvation free energy is given by: 



N 
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(2) 



METHODS 



Model 



The MD simulations and part of the analysis of the 
trajectories were performed with the CHARMM pro- 
gram js^. The oligomeric peptide systems were mod- 
eled by explicitly considering all heavy atoms and the 
hydrogen atoms bound to nitrogen or oxygen atoms 
(PARAM19 potential function ^J). The remain- 
ing hydrogen atoms are considered as part of the car- 
bon atoms to which they are covalently bound (extended 
atom approximation). The effective energy, whose neg- 
ative gradient corresponds to the force used in the dy- 



namics, is 



E{v) = E^acuo{-r) + Gsoiv[-r) 



(1) 



for a molecular system having N heavy atoms with Carte- 
sian coordinates r = (ri, rjv). ^i(r) is the solvent- 
accessible surface computed by an approximate analyti- 
cal expression and using a 1.4 A probe radius. The 
solvation model contains only two a parameters: one for 
carbon and sulfur atoms {ac,s = 0.012 kcal/mol A^), 
and one for nitrogen and oxygen atoms (aAr_o=— 0.060 
kcal/mol A^) [s^- Hence, according to Eq. |21 hydropho- 
bic side chains tend to be buried within the solute 
whereas hydrophilic side chains and the polar groups of 
the backbone prefer to be solvent accessible. Further- 
more, ionic side chains were neutralized |4lj | and a linear 
distance-dependent screening function (e(ry) = 2ry ) was 
used for the electrostatic interactions. The CHARMM 
PARAM19 default cutoffs for long range interactions 
were used, i.e., a shift function [2a| was employed with 
a cutoff at 7.5 A for both the electrostatic and van der 
Waals terms. This cutoff length was chosen to be con- 
sistent with the parameterization of the force-field and 
implicit solvation model. The model is not biased toward 
any particular secondary structure type. In fact, exactly 
the same force field and implicit solvent model have been 
used recently in MD simulations of aggregation [s^, |^ , 
folding of structured peptides (a-helices and /3-sheets) 
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TABLE I: Simulations performed 

Peptide Length T Method IP aggregation lA aggregation 



sequence (m^O (K) events events 

GNNQQNY 10 x 0.5 275 CTMD 6 (19.2) " 

GNNQQNY 5 x 1.0 296 CTMD 3 (14.4) 5 (1.6) 

GNNQQNY 10 x 3.4 330 CTMD 54 (7.6) 43 (1.4) 

GNNQQNY 2 x 1.0 371 CTMD 

GNNQQNY 6 x 2.0 275-400 REMD 14 (60.3) 15 (3.9) 

QQQQQQQ 6 x 2.0 275-400 REMD 27 (54.8) 2 (9.4) 

AAAAAAA 6 x 1.0 275-400 REMD 4 (0.8) 12 (0.9) 

SQNGNQQRG 6 x 2.0 275-400 REMD 1 (1.6) 6 (1.0) 



"The average time (ns) the three peptides remained aggregated 
in IP and lA is given in parentheses 

ranging in size from 15 to 31 residues E3, E3) Q|> and 
small proteins of about 60 residues 14^ . 

REMD simulations 

The basic idea of REMD is to simulate different copies 
{replicas) of the system at the same time but at different 
temperatures values. Each replica evolves independently 
by MD and every tgwap states i,j with neighbor temper- 
atures are swapped (by velocity rescaling) with a proba- 
bility w^j = exp(-A), where A={p,- (3j){Ej - E,), 
f3 = 1/kT and E is the effective energy (potential and 
solvation energy, Eq. A tswap of 10000 MD steps 
(20 ps) was chosen in order to allow the kinetic and po- 
tential energy of the system to relax. High temperature 
simulation segments facilitate the crossing of the energy 
barriers while the low temperature ones explore in de- 
tail energy minima. The result of this swapping between 
different temperatures is that high temperature replicas 
help the low temperature ones to jump across the energy 
barriers of the system. 

In this study six replicas were used with temperatures 
(in K): 275, 296, 319, 344, 371, 400. This range corre- 
sponds to a subset of values used in a previous study of 
reversible peptide folding with the same force-field and 
solvation model The acceptance ratios of exchange 
between neighbor temperatures ranged between 15% and 
24%. Each trajectory has a length of 2 jis for a total of 
12 /is of simulation time (see Tabled). 

Constant temperature MD simulations 

A series of control runs were performed at constant 
temperature: (i) Ten simulations at 330 K (total of 34 
^s) used as a comparison for the aggregation process be- 
tween CTMD and REMD (see Table D); (ii) ten 0.5 ^is 
simulations at 275 K and (iii) five 1 /is simulations at 
296 K to compare CTMD and REMD sampling at phys- 
iologically relevant conditions; (iv) two 1 /xs simulations 
at 371 K to study the system near the condensation tem- 
perature (see below). 



For both REMD and CTMD, Langevin dynamics with 
a friction value of 0.15 ps"^ was used. This friction co- 
efficient is much smaller than the one of water (43 ps^^ 
at 330 K computed as Sirijci/m, where ij is the vis- 
cosity of water at 330 K, and d and m are the effective 
diameter, i.e., 2.8 A , and mass of a water molecule, re- 
spectively) to allow for sufRcient sampling within the fj,s 
time scale of the simulation. The small friction does not 
influence the thermodynamic properties of the system. 

The SHAKE algorithm ^ was used to fix the length 
of the covalent bonds involving hydrogen atoms, which 
allows an integration time step of 2 fs. Furthermore, the 
nonbonded interactions were updated every 10 dynamics 
steps and coordinate frames were saved every 20 ps for 
a total of 5 • 10** conformations//.ts. A 1 run requires 
approximately 2 weeks on a 1.4 GHz Athlon processor 
and the REMD simulations were run in parallel on a 
Linux Beowulf cluster. 

Progress variables 

Aggregation contacts. In-register parallel and an- 
tiparallel aggregation contacts were defined following the 
prescription given in Ref. |30|: a contact was considered 
to be present if the distance between two Ca atoms placed 
on different in-register strands was within 5.5 A. The 
fraction of in-register parallel contacts Qp and in-register 
antiparallel contacts Qa were used to monitor the evolu- 
tion of the aggregation process. In-register parallel and 
antiparallel aggregates, IP and lA respectively, were con- 
sidered formed when Qp and Qa were larger than 0.75 
{Qp,Qa > 11/14) whereas at values smaller than 0.25 
{Qp,Qa < 4/14), the system was considered disordered. 
The aggregation time is defined as the temporal interval 
between the first time point where Qp, Qa < 0.25 and the 
following time point where Qp, Qa > 0.75. 

Radius of gyration. The radius of gyration of the 
oligomeric system Rg was considered to monitor the de- 
gree of condensation and calculated using the minimum 
image convention. Large values of Rg indicate conforma- 
tions with isolated and non interacting peptides {uncon- 
densed phase). Small values of Rg indicate ordered as 
well as disordered aggregated conformations [condensed 
phase). 

Orientational order parameters 

The nematic and polar order parameters, P2 and Pi re- 
spectively, were considered in this study. These order pa- 
rameters represent the first and second rank coefhcients 
of the singlet orientational distribution expanded in a 
Wigner series |4^|H3|i i-^-j a basis set of the Wigner ro- 
tation matrices. The nematic and polar order parameters 
are widely used for studying the properties of anisotropic 
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fluids such as liquid crystals [sil Is^ . Is^ Is^ and are de- 
fined as 

i=l 

and 

_ 1 ^ 

^i = ]^E^-d, (4) 

1=1 

where d (the director) is a unit vector defining the pre- 
ferred direction of alignment, is a suitably defined 
molecular vector, and N is the number of molecules in 
the simulation box, i.e., three peptides in this study. The 
director is defined as the eigenvector of the ordering ma- 
trix [s^ that corresponds to the largest eigenvalue. Here, 
the molecular vectors Zj were defined as unit vectors link- 
ing the peptide's termini (from the N to the C terminus. 
Fig. To optimally select the vectors, other choices 
were investigated: vectors linking the carbonyl C to the 
amide N of each residue ( "amide" vectors) as well as vec- 
tors lying along the carbonyl bonds. Similar results were 
obtained with the three different choices of z^. However, 
due to the atomic connectivity along the backbone the 
"amide" vectors are not fully independent; their orienta- 
tions are strongly correlated and the description of the 
ordered macrostates results less precise. The same is true 
for the "carbonyl" vectors. Hence, vectors linking pep- 
tide's termini were preferred. 




FIG. 1; Pictorial representation of the molecular vectors Zi 
(black arrows) used to compute the order parameters Pi and 
P2. Zi vectors are defined as full length peptide vectors 
(linking the peptide's termini) and allow to clearly discrimi- 
nate between ordered (left, P2 = 0.87) and disordered (right, 
P2 ~ 0.46) conformations of the system. (The pictures were 
drawn using the program PyMOL ||65|)- 

The order parameters (Eq. |31 ^ change value on go- 
ing from one order macrostate to the other and should 
vanish when the transition to a fully isotropic state takes 
place. They describe different orientational properties of 
the system and yield useful and complementary informa- 
tion. The nematic P2 describes the orientational order of 
the system and discriminates between ordered and disor- 
dered conformations. The polar Pi describes the polarity 



of the system, i.e., how much the molecular vectors (z^) 
point in the same direction, and discriminates between 
parallel and antiparallel/mixed ordered aggregates. 

Peptides 

To evaluate the reliability of amyloidogenic propen- 
sity estimations, four oligomeric peptide systems were 
considered in this study: the amyloid-forming heptapep- 
tide GNNQQNY and the soluble nonapeptide SQNGN- 
QQRG both from the yeast prion Sup35 (residues 7-13 
and 17-25 with the Gln/Arg mutation at position 24, 
respectively) l2l|, the amyloidogenic poly(L-glutamine) 
QQQQQQQ |5y| and the non amyloidogenic poly(L- 
alanine) AAAAAAA 57|. To reproduce the experimen- 
tal conditions |21i 1^ , the peptide systems derived 
from the yeast prion Sup35 were modeled without block- 
ing groups, while the Ala and Gin repeats were both 
N-acetylated and C-amidated. 

All simulations were performed with three peptide 
replicas starting from random conformations, positions, 
and orientations. In the initial random positions there 
was no intermolecular contact, i.e., the peptides were sep- 
arated in space. Each system was simulated in a cubic 
box of 75 A per side yielding a sample concentration of 
0.012 M. Since the oligomeric systems present different 
molecular weights, the above reported concentration cor- 
responds to 3.4, 3.9, 5.4 and 3.4 mg/ml for GNNQQNY, 
SQNGNQQRG, QQQQQQQ and AAAAAAA, respec- 
tively. 

Analysis tools 

The aggregation contacts, radius of gyration and order 
parameters analysis was carried out with a GPL licensed 
program |^ developed in house to manipulate and ana- 
lyze MD trajectories. The program is optimized for speed 
and ease of usage so that it allows extensive processing 
of large amounts of data and straightforward addition 
of new analysis tools. Compared to other available pro- 
grams [sil Isgj. the analysis of MD trajectories is much 
faster. 

RESULTS AND DISCUSSION 
REMD diagnostics 

The set of temperatures used in a REMD simulation 
is crucial for a correct and efficient sampling |^ . Since a 
simple "a priori protocol for selecting the optimal tem- 
perature distribution has not been identified (vet), the 
choice often follows empirical considerations U Q, : 
the highest temperature of the set has to be high enough 



5 



i 



10"^ I J — I L — LL — 3 LL^ii — 1 iX 1 — I 

-200 -150 -100 -50 50 

E[kcal/mol] 

FIG. 2: Probability distribution of the effective energy for the 
REMD (solid lines) and the CTMD control simulations (filled 
circles, filled triangles, empty circles and triangles for 275, 
296, 330 and 371 K, respectively). The REMD distributions 
correspond to the following temperatures (from left to right): 
275, 296, 319, 344, 371, and 400 K. The asymmetry of the 
curves and the temperature dependence of the distributions 
indicate the presence of a phase transition around 371 K (see 
text). 
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to overcome energy barriers, while the lowest tempera- 
ture has to allow the exploration of minima. However, 
given a fixed number of replicas the temperature range 
cannot be too wide. Temperature values need to be 
close enough to make the energy histograms overlap (see 
Fig. in order to guarantee a high number of tempera- 
ture swaps during a simulation run. In this study, a set 
of six temperature values ranging from 275 to 400 K has 
been selected (see Methods). The time series of temper- 
ature exchanges for one of the six replicas is shown in 
Fig. O During the simulation, each replica visits all the 
temperatures of the set several times realizing the desired 
free random walk in temperature space 

Symbols in Fig. |21 show the results from CTMD simu- 
lations carried out at 275 (filled circles), 296 (filled trian- 
gles), 330 (empty circles) and 371 K (empty triangles). 
At 330 K, the CTMD effective energy distribution is 
located between the REMD distributions extracted at 
319 and 344 K and shows a consistent functional pro- 
file. At 371 K, CTMD and REMD effective energy dis- 
tributions overlap. Therefore, the energetic properties of 
an aggregating system sampled by a REMD simulation 
at medium and high temperatures correspond to those 
observed in CTMD simulations. However, approaching 
the physiologically relevant conditions the CTMD dis- 
tributions tend to shift toward less favorable energies 
(Fig 121 filled symbols). CTMD at low temperature can 
get trapped in local energy minima and REMD is supe- 
rior in sampling conformational space 0, . 

The time series of the fraction of in-register paral- 
lel contacts (Qp) and in-register antiparallel contacts 



FIG. 3; Time series of (from top to bottom) the temperature 
T, the fraction of in-register parallel contacts Qp, and the 
fraction of in-register antiparallel contacts Qa for a REMD 
replica. Along the trajectory, replicas realize the desired free 
random walk in temperature space (top) so that an efficient 
sampling of the ordered aggregates is allowed (peaks in Qp 
and Qa plots). Horizontal lines in the time series of the frac- 
tion of aggregation contacts indicate the upper/lower thresh- 
olds used to define the ordered aggregation/disaggregation 
events. 

{Qa) have been monitored along the REMD trajectories 
(Fig. 13). A total of 14 IP and 15 lA aggregation events 
have been observed along the total simulation time of 
12 /is (see Tabled). The average aggregation time (see 
Methods) was 0.74 fj,s for IP and 0.75 fis for lA arrange- 
ments. The average aggregation time determined from 
the REMD simulation is similar to the values obtained 
from 34 /xs CTMD simulations at 330 K. It is worth not- 
ing that in a preliminary REMD run with higher tem- 
peratures values (6 x l^s, 319-465 K; data not shown) 
only 3 IP and 4 lA aggregation events were sampled. 
The temperature range is crucial in REMD and it has 
to be carefully chosen in order to smed up the confor- 
mational search of relevant states |6(l |. i.e., the ordered 
states when studying aggregation. To bias the search to- 
ward conditions where ordered states are more probable, 
the temperature was set to lower values (275 to 400 K as 
mentioned above) and the sampling of aggregation events 
turned out substantially improved. 

Fig. 21 shows the projections of the free energy surface 
along Qp and Qa for both REMD and CTMD trajecto- 
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FIG. 4: Free-energy projections along the fraction of in- 
register parallel contacts Qp (left) and in-register antiparal- 
lel contacts Qa (right). Conformations with zero in-register 
contacts were chosen as reference states. AG'(„_q) was com- 
puted as —kBT\n{N„/No), where N„ indicates the number 
of conformations with n contacts and fc_g is the Boltzmann 
constant. REMD data are shown in solid lines for all the 
temperature values except for 319 and 344 K which are in 
dashed lines. CTMD data are shown with symbols (filled cir- 
cles, filled triangles, empty circles and triangles for 275, 296, 
330 and 371 K, respectively). 

ries. The profiles indicate that the structural properties 
of the aggregating system sampled by a REMD simu- 
lation correspond to those observed in CTMD simula- 
tions only at high and medium temperatures. At 371 K, 
CTMD and REMD free energy projections overlap. At 
330 K, the CTMD free energy profiles (empty circles) 
are correctly placed between REMD projections at 319 
and 344 K (dashed lines) and show patterns character- 
ized by a well defined local minimum at Qp > 0.7 and 
a monotonic uphill trend along Qa, fully consistent with 
the profiles extracted from the REMD simulation. How- 
ever, at low temperature (275 and 296 K) the free energy 
profiles extracted from CTMD and REMD trajectories 
are not consistent any more and the most "relevant" con- 
formations, which correspond to in-register parallel and 
antiparallel arrangements {Qp,Qa > 0.7), are not cor- 
rectly sampled by CTMD (Fig El filled symbols). 



Temperature dependence of ordered amyloid 
peptide aggregation 

Since the energetic and structural properties of the sys- 
tem are not artificially altered (see previous subsection), 
the REMD approach allows to evaluate thermodynamic 
quantities as a function of temperature in the chosen 
range 6]. From the REMD simulation performed for 
this study, the properties of interest have been extracted 
at any temperature of the set and the aggregation of 



the amyloid-forming peptide GNNQQNY has been mon- 
itored in temperature space (275-400 K). This analysis 
gives interesting insights into the amyloid aggregation 
process. 

The effective energy histograms shown in Fig. |21 are 
not symmetrically distributed around their mean value 
and their shape varies with temperature. The distribu- 
tions, in fact, broaden toward higher energy values at 
low temperature (275-344 K) and toward lower energy 
values at high temperature (371-400 K). Moreover, by 
increasing temperature they progressively become lower 
and broader till the value of 371 K is reached. Mitsutake 
et al. have interpreted such a behavior as the evidence 
of a phase transition [6]| . To characterize the transi- 
tion, the radius of gyration Rg of the oligomeric system 
was considered and free energy projections along Rg were 
plotted (see Fig. [SJ. Conformations of the system pro- 
ducing non interacting peptides, namely conformations 
where all inter-peptide atomic distances are larger than 
the long-range interactions cutoffs (7.5 A in this case), 
were used to determine i?g*", i.e., the lowest detected ra- 
dius of gyration for isolated peptides (see Fig. [SJ. The 
existence of two macrostates in equilibrium has been re- 
vealed: the first, named uncondensed state, includes high 
energy conformations with one or more isolated peptides 
{Rg > Rg''' ); the second, named condensed state, con- 
sists of low energy conformations with aggregated pep- 
tides {Rg < Rg^)- For entropic reasons, the uncondensed 
state is preferred at high temperature. By cooling down, 
the condensed state is increasingly stabilized and around 
371 K the fiuctuations of Rg show a well defined peak 
highlighting the presence of the condensation transition 
(see Fig. |5l . The equilibrium between the condensed and 
the uncondensed macrostates is clearly concentration de- 
pendent. If the concentration of amyloid- forming units 
increases, the equilibrium is moved toward the condensed 
state and the aggregation process is favored. 

The free energy profiles along Qp and Qa at various 
temperatures help in understanding how the nucleation 
process evolves upon peptides condensation. At val- 
ues of 400, 371 and 344 K both projections show steep 
uphill patterns with a single free energy minimum at 
Qp ^ Qa ^ (see Fig. 0}. This means that upon con- 
densation the peptides are still more likely to form dis- 
ordered aggregates characterized by non-specific interac- 
tions than amyloid-forming nuclei. In this range of tem- 
peratures, the enthalpic contribution due to in-register 
backbone or side-chain interactions does not dominate 
the entropic one and the growth of ordered nuclei is for- 
bidden. However, when the temperature decreases the 
entropic contribution becomes less important and or- 
dered in-register aggregates start forming. As shown in 
Fig. El in fact, below 330 K two and one additional free 
energy minima appear in the projection along Qp and 
Qa, respectively. The observed minima correspond to 
in-register parallel {Qp > 0.7) and in-register mixed or 
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out-of-register (0.3 < Qp < 0.7 and 0.4 < Qa < 0.7) ar- 
rangements and strongly suggest that the three-peptide 
system moves toward a higher degree of order when ap- 
proaching the physiologicaUy relevant conditions. 

The simulation results indicate that in the early steps 
of amyloid aggregation a condensation stage toward dis- 
ordered aggregates precedes the nucleation process and 
the disorder-order transition, in agreement with experi- 
mental evidence 1621. 
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FIG. 5: (Top) Free-energy projections along the radius of gy- 
ration of the oligomeric system Rg computed from REMD 
trajectories. Solid lines correspond to temperature values be- 
low the condensation temperature (275-344 K); dashed lines 
correspond to temperature values above the transition tem- 
perature (371 and 400 K). The lowest radius of gyration for 
the uncondensed state is shown as a vertical line {Rg'~^ — 11.9 
A ). (Bottom) Temperature dependence of the average radius 
of gyration (Rg) (filled circles) and its fluctuations a (empty 
circles). The behavior of {Rg} and a indicates the presence 
of a phase transition around 371 K between a condensed (low 
T) and an uncondensed phase (high T). Fluctuations of the 
radius of gyration a are computed as (-Rg) — (Rg)'^- Data at 
431 and 465 K were obtained from a preliminary REMD run 
carried out in a higher temperature range (6 x Ifis, 319, 344, 
371, 400, 431 and 465 K). 



Disorder-order transition 

In the early steps of aggregation, amyloidogenic pe^ 
tides assemble into highly ordered /3-sheet structures 
l30| . During the assembly, the peptides tend to align 
adopting an extended /3-strand conformation and a re- 
markable change in the local orientational order occurs. 
The aggregation of amyloid-forming peptides may then 
be interpreted as an order transition and orientational 
order parameters are suitable to monitor the time evolu- 
tion of the process. Two orientational order parameters 
were employed and free energy projections are shown in 
Fig. IHI Along P2, the free energy profiles show a first 
broad minimum at P2 ~ 0.5 for any temperature of the 
set and a second narrower one at P2 « 0.9 for T values 
below 330 K. The first corresponds to a large free en- 
ergy basin where orientational order is absent while the 
second corresponds to a smaller and well defined basin 
with a high orientational degree of order. Although the 
order parameters should vanish when order is absent. 
Fig. IHI shows that this is not the case when the number 
of vectors is small. Since only three peptides were sim- 
ulated, a "background" order was always detected and 
the free energy minimum describing the disordered state 
is placed at P2 « 0.5 which is consistent with the value 
of y/Sl/AOir N expected for a completely randomly ori- 
ented array of N molecules The order parameter P2 
shows the existence of two macrostates in equilibrium: 
the disordered state with a high entropy content, which 
corresponds to the global minimum of the free energy 
surface at high temperature, and the ordered state which 
becomes the global free energy minimum at low tempera- 
ture. Interestingly, the free energy profiles along Qp and 
Qa do not lead to the same conclusion and the observed 
in-register arrangements correspond to local minima of 
the free energy surface (see Fig.^J. 

Along Pi, two narrow and well distinct minima cor- 
responding to ordered macrostates at different polarity 
appear on the free energy projections (Fig.lHJ. The first, 
displayed at Pi « 0.35, describes a free energy basin with 
a high-order and low-polarity content. Conversely, the 
second, displayed at Pi « 0.95, corresponds to a basin 
with a high-order and high-polarity content. The order 
parameter Pi discriminates between parallel and antipar- 
allel/mixed ordered conformations and provides comple- 
mentary information since it allows to further character- 
ize the ordered state. 

Symbols in Fig. El show the free energy projections 
along the order parameters from CTMD simulations. 
Once again, the comparison with REMD profiles indi- 
cates that isothermal MD (filled symbols) does not sam- 
ple the ordered aggregates with their correct statistical 
weight close to the physiological temperature range. 

The REMD free energy profiles along Pi show that at 
low temperature (275 and 296 K) both polar macrostates 
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are highly populated. In the investigated temperature 
range, the system does not show an overall polar degree 
and frequent jumps between ordered states characterized 
by different polarity are observed. This suggests that be- 
low the order transition the equilibrium between polar 
macrostates might help amyloidogenic systems overcom- 
ing the entropy loss occurring during nucleation. In other 
words the growth of amyloid-forming nuclei might have 
an entropically favorable component due to the multiple 
ordered macrostates. 



Sequence dependence of amyloidogenic propensity 

Free energy projections along the nematic order pa- 
rameter P2 show how the equilibrium between the or- 
dered and disordered state changes in temperature space 
(Fig. EJ. Upon cooling, the statistical weight of the or- 
dered state increases and the mean of the P2 distribution 
moves toward higher values. The value of (^2), where 
(• • • ) indicates the average over the canonical ensemble, 
is then related to the thermodynamic stability of the 
ordered state and could be used to measure the amy- 
loidogenic propensity of the system. (P2) values com- 
puted at different temperatures from REMD trajectories 
of the amyloid-forming peptide GNNQQNY are shown 
in Fig. [71 with filled circles. At high temperature, the 
(P2) values are close to 0.5 because no orientational or- 
der is present, and the system does not show amyloido- 
genicity. By decreasing temperature, the amyloidogenic 
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P2 Pi 

FIG. 6: Free-energy projections along the nematic (P2, left) 
and the polar (Pi, right) order parameters. REMD data are 
shown in solid lines for all the temperature values except for 
319 and 344 K which are in dashed lines. CTMD data are 
shown with symbols (filled circles, filled triangles, empty cir- 
cles and triangles for 275, 296, 330 and 371 K, respectively). 
Schematic representations of the aggregates (black arrows) 
are depicted to show that order parameters yield complemen- 
tary information: P2 discriminates between ordered and disor- 
dered conformations while Pi discriminates between parallel 
and antiparallel/mixed ordered aggregates. 



0.85 




T(K) 

FIG. 7: Temperature dependence of the nematic order pa- 
rameter {P2) averaged over the canonical ensembles sampled 
by REMD for four oligomeric peptide systems. (P2) estimates 
the amyloidogenic propensity of peptide systems and discrim- 
inates between amyloidogenic (GNNQQNY and QQQQQQQ) 
and non amyloidogenic (SQNGNQQRG and AAAAAAA) se- 
quences in agreement with experimental data plLlSdls^ . 

propensity grows and becomes increasingly larger until 
the order transition is completed. At physiologically rel- 
evant conditions, (P2) « 0.65 and the system is hi ghly 
amyloidogenic in agreement with experimental data [21j . 

Since the orientational order parameters do not de- 
pend on the peptide sequence and length, the reliability 
of the predictions could be further tested in sequence 
space. The REMD protocol was then applied to three 
additional oligomeric peptide systems (see "Methods") 
and (P2) values were evaluated to measure and compare 
amyloidogenic propensities. The testing set comprises 
a nonapeptide from the yeast prion Sup35 (SQNGN- 
QQRG) experimentally studied by Balbirnie et al. 
and two heptapeptides (QQQQQQQ and AAAAAAA). 
Glutamine and alanine homopolymers flanked by basic 
residues to iinprove solubility have been investigated by 
Perutz et al. IsilsTf. 

Experimentally, the nonapeptide SQNGNQQRG 
shows solubility in vivo and in vitro and no formation 
of amyloid fibrils |2]| . In agreement with these findings, 
(P2) is smaller than 0.55 in the whole temperature range 
(Fig. 13 empty squares) and the system is considered as 
non-amyloidogenic. The number of aggregation events 
and the average life time of aggregation extracted from 
REMD trajectories and reported in TableQ] Remarkably, 
these quantities show that non-amyloidogenic sequences, 
i.e, SQNGNQQRG and AAAAAAA, do transiently as- 
semble in a /3-sheet conformation but still remain soluble 
because their ordered aggregates do not correspond to 
well defined free energy minima. 

CD spectra, electron micrographs and X-ray diffrac- 
tion photographs showed that poly(L-glutamine) pep- 
tides aggregate in solution at both pH 7.0 and 3.0 form- 
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FIG. 8: (Top) Snapshots of ordered aggregates of three 
(thick sticks) and six (thin sticks) amyloidogenic SYVIIE pep- 
tides Isl extracted from CTMD simulations at 330 K. The 
simulations were performed at a sample concentration of 5 
mg/ml. The overall conformation and twist of the three- 
stranded and six-stranded parallel /3-sheets are indistinguish- 
able. (Bottom) The six-stranded /3-sheet upon 90° rotation 
to better visualize the twist. (The pictures were drawn using 
the program PyMOL 



CONCLUSIONS 

The present study shows that atomistic REMD sim- 
ulations with implicit solvent allow to sample the early 
steps of ordered aggregation of amyloidogenic peptides 
at physiologically relevant temperatures. The free en- 
ergy profiles projected along structural and orientational 
progress variables are essentially the same in REMD and 
CTMD. The discrepancies at temperature values below 
330 K are due to the limitations in sampling in CTMD 
simulations which indicates that REMD is a more effi- 
cient approach in the physiological range. 

The early steps of amyloidosis can be interpreted as 
a condensation followed by an order transition. There- 
fore, the REMD simulation results were analyzed with 
two order parameters originally introduced to study liq- 
uid crystals. Interestingly, the nematic order parameter 
averaged over a canonical ensemble is able to discrimi- 
nate amyloidogenic from soluble peptides in agreement 
with experimental data. 

Although the present study was performed with three 
peptides for reasons of computational efficiency, the 
description of the ordered aggregates is likely to be 
independent of the size of system, i.e., the number of 
simulated peptide replicas. Very recent MD simulations 
of the amyloidogenic SYVIIE peptide (Cecchini et ai, 
unpublished results), which has been experimentally 
investigated by de la Paz and Serrano |64l| . have shown 
ordered aggregates of six peptides. Interestingly, the 
parallel /3-sheet consisting of six peptides has the same 
overall conformation and twist as the three-peptide 
aggregate (Fig. |SJ). 



ing tightly linked /3-sheets structures (5^. In particular, 
the X-ray diffraction picture exhibits a fiber diagram of 
the cross-/3 type distinctive of amyloid fibrils. On the 
other hand, poly(L-alanine) doesn't display amyloido- 
genicity and CD spectra showed a-helical structures at all 
pHs [53 ■ Again, the (P2) patterns shown in Fig. [7| (filled 
squares and empty circles) are consistent with experi- 
mental findings and correctly indicate amyloidogenicity 
only for QQQQQQQ. 

Interestingly, Fig. [7| allows also to compare between 
amyloidogenic sequences. In fact, according to the (P2) 
patterns the glutamine repeat is more amyloidogenic 
than GNNQQNY at physiologically relevant conditions. 
To our knowledge, no experimental data is available to 
verify this finding. Testing of this prediction is a chal- 
lenge for experimentalists. 
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